Variational quantum evolution equation solver

Variational quantum algorithms offer a promising new paradigm for solving partial differential equations on near-term quantum computers. Here, we propose a variational quantum algorithm for solving a general evolution equation through implicit time-stepping of the Laplacian operator. The use of encoded source states informed by preceding solution vectors results in faster convergence compared to random re-initialization. Through statevector simulations of the heat equation, we demonstrate how the time complexity of our algorithm scales with the Ansatz volume for gradient estimation and how the time-to-solution scales with the diffusion parameter. Our proposed algorithm extends economically to higher-order time-stepping schemes, such as the Crank–Nicolson method. We present a semi-implicit scheme for solving systems of evolution equations with non-linear terms, such as the reaction–diffusion and the incompressible Navier–Stokes equations, and demonstrate its validity by proof-of-concept results.

where n x , n y and n t are prescribed positive integers, such that x i = x L + i · x , y j = y L + j · y , t k = k · t , �x = L x /n x , �y = L y /n y and �t = T/n t , where L x = x R − x L and L y = y R − y L . The discrete domain grid is denoted by � d = {(x i , y j ) : n x ∈ {0, 1, . . . , n x }, n y ∈ {0, 1, . . . , n y }} and boundary grid by Ŵ d .
The finite difference (FD) approximation for the second-order spatial derivative (5-point) of the Laplacian operator taken at t = t k is where δ x := D�t/�x 2 and δ y := D�t/�y 2 are diffusion parameters.
Using first-order FD for temporal derivative (u k+1 ij − u k ij )/�t weighted by ϑ ∈ [0, 1] , the evolution Eq. (1) can be expressed in vector shorthand as where I is the identity matrix of the same size, u k = u k ij 0≤i≤n x ,0≤j≤n y and f k = f k ij 0≤i≤n x ,0≤j≤n y . Depending on the choice of parameter ϑ , actual time-stepping may follow an explicit (forward Euler) method ( ϑ = 0 ), an implicit (backward Euler) method ( ϑ = 1 ), a semi-implicit (Crank-Nicolson) method ( ϑ = 1/2 ) or a variable-ϑ method 24 . The explicit method ( ϑ = 0 ) is efficient for each time-step but is only stable if it satisfies the stability condition v ≤ 1/2 . The implicit (backward Euler) method ( ϑ = 1 ) is unconditionally stable and first-order accurate in time ( ε ∼ �t ), which reads The semi-implicit Crank-Nicolson (CN) method ( ϑ = 1/2 ) is popular as it is not only stable, but also secondorder accurate in both space and time ( ε ∼ �t 2 ), which reads where f k+1/2 = (f k+1 + f k )/2 . However, the CN method may introduce spurious oscillations to the numerical solution for non-smooth data unless the algorithm parameters satisfy the maximum principle 25 . Variational quantum solver. Here, we explore a variational quantum approach towards the solution of the evolution equation (1). In addition to potential quantum speedup, a variational quantum algorithm could also benefit from data compression, where a matrix of dimension N can be expressed by a quantum system x ij , t k := x i , y j , t k , i = 0, 1, . . . , n x ; j = 0, 1, . . . , n y ; k = 0, 1, . . . , n t , (6) Au i,j = −δ x (u i−1,j − 2u i,j + u i+1,j ) − δ y (u i,j−1 − 2u i,j + u i,j+1 ), www.nature.com/scientificreports/ with only log 2 N qubits, where N is the size of the problem. Consider the Poisson equation, which is a timeindependent form of Eq. (1), expressed as The Laplacian operator ∇ 2 in one dimension can be discretized using the finite difference method in the x direction into an N N coefficient matrix A x,β as where β ∈ {D, N} refers to either the Dirichlet (D) or Neumann (N) boundary condition, and α D = 1 and α N = 0 . This extends naturally to higher dimensions, for instance A y,β in the y direction.
A variational quantum solution is to prepare a state |u� such that A|u� is proportional to a state |b� in a way that satisfies Eq. (10). To do that, a canonical approach 10-12 is to first decompose the matrix A over the Pauli basis P n = {P 1 ⊗ · · · ⊗ P n : ∀i, P i ∈ {I, X, Y , Z}} (where X = |1��0| + |0��1| , Y = i|1��0| − i|0��1| , and Z = |0��0| − |1��1| are the Pauli matrices and I = |0��0| + |1��1| is the identity matrix) as where c P = tr (PA)/2 n are the coefficients of A in the Pauli basis. Using simple operators σ + = |0��1| , σ − = |1��0| , the number of terms in the decomposition can be reduced to 2 log N + 1 17 . A more efficient approach, however, is to express A as a linear combination of unitary transformations of simple Hamiltonians 18 . Accordingly, the decomposition of A in one dimension can be written as 19 where I 0 = |0��0| and β ∈ {D, N} , as before except here, a D = 0 and a N = 1 . Here, S is the n-qubit cyclic shift operator defined as The expectation values of a Hamiltonian H including the shift operator S are evaluated by applying the unitary shift operator to the quantum state 18 , where |φ� is an arbitrary n-qubit state and |φ ′ � = S|φ� . Note that Eq. (13) can be re-written as Since expectation values of the identity operator are equal to 1, i.e. �φ|I ⊗n |φ� = �φ ′ |I ⊗n |φ ′ � = 1 , evaluating the expectation value of the operator A x,β requires only the evaluation of expectation values of the simple Hamiltonians H 1−4 ( H 1−3 for Dirichlet boundary condition). The required number of quantum circuits is therefore limited to a constant O(n 0 ) 18 . Similar decomposition expressions apply to problems of higher dimensions, including A y,β in the y direction 19 .
Once the matrix A is decomposed, a parameterized quantum state |ψ(θ)� is prepared using an Ansatz represented by a sequence of quantum gates U(θ) parameterized by θ applied to a basis state |0� ⊗n , such that |ψ(θ)� = U(θ)|0� ⊗n . Here, we use a hardware-efficient Ansatz consisting of multiple layers of R Y gates across n qubits entangled by controlled-X gates (see Fig. 1). For the source term f in (10), a quantum state |b� is prepared by encoding a real vector with the unitary U b , such that |b� = U b |0� ⊗n . Depending on the actual input, conventional amplitude encoding methods 26,27 may introduce a global phase that must be corrected by its complex argument for computing in the real space.
Using classical optimization tools, the cost function (17) is minimized with θ updated iteratively until convergence is reached. The optimization process follows either a gradient-based or gradient-free approach, depending on how the gradient of the cost function is evaluated. A gradient-free optimizer is guided by an estimate of the inverse Hessian matrix, whereas a gradient-based optimizer by the partial derivative of the cost function E with respect to parameters θ , i.e. ∂E/∂θ , which can be evaluated by a quantum computer (for details, see 18,19 ). Regardless of the choice of gradient optimizer used, the optimization routine halts when the cost error falls under a convergence threshold ( ǫ < ǫ tol ) whence the parameters are at optimum θ = θ opt . The converged solution vector |u� = r opt |ψ(θ opt � satisfies 10 where r opt = r(θ opt ) is the norm of the solution to the Poisson equation (10).
In this study, we propose to solve the evolution equation (1) through successive time-stepping of the quasisteady Poisson equation using a variational quantum algorithm. Using a parameter set θ k obtained at time-step k, we encode a normalized source state |b k � := |b�/ √ �b|b� from |b(θ k )� and seek an implicit solution to where θ k+1 = θ opt (t k+1 ) is the parameter set and r k+1 = r opt (θ k+1 ) is the norm at next time-step k + 1 . This process is then iterated in time up to n t number of time-steps as desired (see Algorithm 1).
. . . R y (θ n n ) Figure 1. Schematic of the quantum circuit hardware-efficient Ansatz used in this study. www.nature.com/scientificreports/ The cost function E(θ k ) may be evaluated on a quantum computer by computing each of the inner products in the expression in Eq. (17) separately. Using the decomposition provided in Eq. (16), these inner products may be expressed in terms of expectation values �ϕ|O i |ϕ� for preparable states |ϕ� and simple Hermitian operators O i . Each expectation value �ϕ|O i |ϕ� is evaluated on a quantum computer by preparing the state |ϕ� using the quantum circuits described above and then measuring the operator O i in the state |ϕ� 18 .
In this study, the variational quantum algorithm is implemented in Pennylane (Xanadu) 28 using a statevector simulator with the Qulacs 29 plugin as a backend for quantum simulations, and the L-BFGS-B optimizer for parametric updates. Amplitude encoding is carried out via the standard Mortonnen state preparation template 30 with custom global phase correction. For hardware emulation via the QASM simulator (Qiskit), we refer the reader to the excellent cost-sampling analysis of Sato et al. 18 .

Applications to the heat/diffusion equation
Consider the following one-dimensional heat or diffusion equation without a source term Dirichlet conditions are applied on the boundaries of a 1D domain To solve Eq. (22), the variational quantum evolution algorithm (Algorithm 1) can be employed with a suitable time-stepping scheme (7). For the implicit Euler (IE) method (8), the matrix A and source state |b(θ k )� can be decomposed into For the Crank-Nicolson (CN) method (9), it follows that the Au k term, which carries a small but non-trivial evaluation cost, can be eliminated using the source state of the previous time-step k − 1 , leading to Here, the presence of a k − 1 term in b k−1 is not unexpected due to temporal finite differencing at second-order accuracy.
For a space-time domain × J ∈ [0, 1] × [0, 1] , let the number of time-steps be n t = 20 and the spatial intervals be n x = 2 n + 1 , where n is the number of qubits, and δ x = 1 is the diffusion parameter. We employ the Dirichlet boundary condition with boundary values (g L , g R ) = (1, 0) and initial values u 0 = 0 . With initial random parameters θ 0 ∈ [0, 2π] , we run a limited-memory Broyden-Fletcher-Goldfarb-Shanno boxed (L-BFGS-B) optimizer [31][32][33][34] to optimize θ with absolute and gradient tolerances set at 10 −8 .  Figure 2b shows how the cost function E depends on the number of optimization steps for n-l of 3-3 and 4-4 (10 sampled runs each). Each distinct step in E represents sequential optimization from solution |ψ(θ k )� at time-step k towards the solution |ψ(θ k+1 )� at k + 1 . For small time-step t , θ k provides a good initial parameter set for solving optimization step k + 1 . If the Ansatz parameters were re-initialized randomly θ k ∈ [0, 2π] before each time-step, significantly more optimization steps would be required on average for convergence for each run (see Fig. 2b, inset).

Time complexity.
Here we briefly examine the time complexity of the quantum algorithm excluding the classical computing components. Following the analysis of the variational Poisson solver 18 , the time complexity of the proposed variational evolution equation solver per time-step reads where the terms within the inner parentheses indicate the time complexity of the state preparation scaling as O(l + e + n 2 ) , which consists of the Ansatz depth l, the encoding depth e = O(n 2 ) 35 , the depth of the circuit  where T eval is the sum of function evaluations required for a run with n t time-steps. Using a gradient-based optimizer, the time complexity for gradient estimation via quantum computing would scale as the Ansatz volume O(nl) representing the number of quantum circuits required for parameter shifting. Otherwise, with a gradientfree optimizer, the time complexity simply contributes towards T eval as additional function evaluations required to evaluate the Hessian for gradient descent.
To see if the time complexity for gradient-free optimization scales as O(nl) , we plot the time-averaged number of function evaluations T eval against the number of parameters nl (Fig. 3a). Indeed, we found that T eval scales reasonably with nl (see trendline of slope 1), despite apparent tapering at higher l. Figure 3b shows that the timeaveraged trace error ε tr decreases with circuit depth l, even for over-parameterized quantum circuits where the number of layers exceeds the minimum required for convergence, l min := 2 n /n 36 . For low grid resolution n = 3 , the trace error is limited to a minimum of ∼ 10 −4 . Since the time complexity for solving the Poisson equation classically is O (N log 2 N) , where N = 2 n , quantum advantage could be realized with the proposed algorithm with linear time scaling by n t at sub-exponential time complexity 18 . www.nature.com/scientificreports/ For deep and wide quantum circuits, the increase in optimization time is exacerbated by the presence of barren plateaus, or vanishingly small gradients in the energy landscape, where re-initialization can leave one trapped at a position far removed from the minimum [37][38][39] . Conversely, short time-steps lead to efficient solution trajectories that remain close to the local cost minima, leading to significant reduction in optimization times. To verify this, we conduct numerical simulations varying the diffusion parameter δ with l = n up to time T = 1 . Figure 3c shows that the number of iterations, or required optimization steps, per time-step increases linearly with δ . Close inspection of the time-averaged trace distance shows bimodal distributions at higher δ , which separates success and failure during convergence towards the global minimum (see Fig. 3d, dotted lines), resembling local minima traps due to poor optimization or expressivity of Ansätze 19,40 . Discretization error. Time evolution can be at a higher order, specifically for the Crank-Nicolson method.
The problem statement is identical to the previous one, except with Dirichlet boundary values (g L , g R ) = (0, 0) and the initial condition u 0 = sin(πx/L x ), where we use L x = 1 as the spatial length of the domain. This admits an exact analytical solution, Figure 4 compares variational quantum and exact solutions using implicit Euler and Crank-Nicolson (CN) schemes. The discretization error for the higher-order CN scheme is reduced significantly, especially at lower grid resolution ( n = 3 ). Although the complexity costs for both methods (23) and (24) are similar, note however that the CN method may introduce spurious oscillations for non-smooth data 24 , an issue which may be exacerbated by quantum noise 41 .
Under the implicit Euler scheme (8), the matrix A and source state |b k � can be decomposed into Dirichlet conditions are applied on the boundaries, where u(x L,R , y, t) = g x L,R (y, t) and u(x, y L,R , t) = g y L,R (x, t) . Let the number of spatial grid intervals be n x = 2 m x + 1 and n y = 2 m y + 1 , where m x is the number of qubits allocated to the x grid, m y to the y grid, and n = m x + m y is the total number of qubits. Accordingly, A is decomposed in terms of simple Hamiltonians in x and y as where H 1−4 are simple Hamiltonians to be evaluated ( H 1−3 for Dirichlet boundary condition). Figure 5 shows solution snapshots to a 2D heat conduction or diffusion problem taken at time T = 1 with Dirichlet boundary values (g x L , g x R ) = (0, 0) and (g y L , g y R ) = (1, 0) , initial values u 0 = 0 , n t = 20 and diffusion parameters δ x = δ y = 1 . Results obtained from variational quantum solver agree with classical solutions with time-averaged trace errors of up to 10 −2 .

Applications to the reaction-diffusion equations
Here, we extend applications of our variational quantum solver to evolution equations with non-trivial source terms. Consider a two-component homogeneous reaction-diffusion system of equations x,D + δ y u k+1 y,D . Implementation. The semi-implicit variational quantum solver solves for the Laplacian for each component using a quantum computer and the solution vectors are explicitly coupled through source terms prior to re-encoding in preparation for the next time-step (see Algorithm 2).  Figure 6a shows how an initial mid-pulse can spontaneously and periodically split in space and time, a phenomenon captured using variational quantum diffusion reaction solver (see Algorithm 2) even on relatively low spatial resolutions.
1D Brusselator model. So far, we have been looking at only Dirichlet boundary conditions. Here, we demonstrate a test example for Neumann boundary conditions in a diffusion-reaction model, namely, the Brusselator model 48 Figure 6b shows how a chemical pulse can be spontaneously created, which continually travels leftwards in time, creating traveling waves that appear as striped patterns in time despite low spatial resolutions.

Applications to the Navier-Stokes equations
The Navier-Stokes equations are a set of non-linear partial differential equations that describes the motion of fluids across continuum length scales. There are several studies aimed at applying quantum algorithms to computational fluid dynamics (see review 49 ), ranging from reduction of partial differential equations to ordinary differential equations 50 and quantum solutions of sub-steps of the classical algorithm 51,52 to the quantum Lattice Boltzmann scheme 53 .
Here, we look into the potential use of variational quantum algorithms to evolve the fluid momentum equations in time. Consider the incompressible Navier-Stokes equations in non-dimensional form www.nature.com/scientificreports/ where u is the velocity vector and p is the fluid pressure. The ratio Re = U c L c /ν is the Reynolds number, where U c is the characteristic flow velocity across a characteristic length-scale L c and ν is the fluid kinematic viscosity. Unlike other temporal evolution equations, the incompressible Navier-Stokes equations cannot be timemarched directly as the resultant velocities do not satisfy the continuity constraint (45), and hence are not divergence-free. To resolve this, the projection method 54 , also known as the predictor-corrector or fractional step method, separates the solution time-step into velocity and pressure sub-steps, also known as the predictor and corrector steps.
Projection method. Predictor step. The predictor step first approximates an intermediate velocity u * by solving the fluid momentum equation (44) in the absence of pressure, i.e. the Burgers' equations 55 , of the form Through a semi-implicit scheme, the viscous terms are handled implicitly using the variational quantum evolution equation solver and the non-linear inertial terms explicitly as source terms using classical computation. For quantum algorithms for non-linear problems, the reader is referred to separate works on quantum ordinary Implementation. Overall, the variational quantum solver for Navier-Stokes equations using the projection method (see Algorithm 3) involves two sequential steps, the first requiring a number of Algorithm 1 iterations equal to the number of velocity components, and the second for the pressure Poisson step. For two-dimensional flows, the number of velocity components to be solved can be effectively reduced by one through the vorticity stream-function formulation (Appendix B). In computational fluid dynamics, these implicit systems of linear equations are often the most computationally expensive parts to solve in classical algorithms, providing incentives for potential speedup via quantum computing 49,52 .  57 . Consider a two-dimensional square domain � = (0, L) × (0, L) ⊂ R 2 with only one wall sliding tangentially at a constant velocity. For simplicity, we employ a fixed collocated grid, instead of a staggered grid which helps avoid spurious pressure oscillations but at the cost of increased mesh and discretization complexity.
No-slip boundary conditions apply on all walls, so that zero velocity applies on all wall boundaries except one moving at u(x, 0) = 1. Figure 7a shows a snapshot of a test case conducted on a 2 n = 8 × 8 grid at t = 0.5 up to T = 5 , with the central vortex shown by normalized velocity quivers in white. In terms of time complexity, we note that the pressure correction step requires a greater number of function evaluations for convergence compared to an implicit velocity step (Fig. 7b). This is due to the additional quantum circuits for evaluating the H 4 Hamiltonians (32,  www.nature.com/scientificreports/ 33) for Neumann boundary conditions and one for specifying the reference pressure (50), leading to a total of 9 evaluation terms compared to 6, a ratio which corroborates with the apparent ∼ 50% increase in function evaluations shown in Fig. 7b.
While not directly comparable to classical computational fluid dynamics in numerical accuracy, this exercise, nevertheless, roadmaps potential applications of the variational quantum method towards more complicated flow problems 52 .

Conclusion
In this study, we proposed a variational quantum solver for evolution equations which include a Laplacian operator to be solved implicitly. For short time-steps t , the use of initial parameter sets encoded from prior solution vectors results in faster convergence compared to random re-initialization. The overall time complexity scales with the Ansatz volume O(nl) for gradient estimation and with the number of time-steps O(n t ) for temporal discretization. Our proposed algorithm extends naturally to higher-order time-stepping and higher dimensions. For evolution equations with non-trivial source terms, the semi-implicit scheme can be applied, where non-linear source terms are handled explicitly. Using statevector simulations, we demonstrated that variational quantum algorithms can be useful in solving popular partial differential equations, including the reaction-diffusion and the incompressible Navier-Stokes equations. Together, our proposed algorithm extends the use of quantum Poisson solvers to solve time-dependent problems with reduced time complexity from variational quantum algorithms over classical computation.
The present work aims at bridging the gap between variational quantum algorithms and practical applications. Our work has assumed that the state preparations, unitary transformations and measurements are implemented perfectly, and does not consider the effects of quantum noise from actual hardware or any potential amplification from iterative time-stepping. In our implementation, we only considered the hardware-efficient Ansatz with R y rotation gates and controlled-NOT entanglers, and thus leave open the question about the performance of other Ansätze. Future work can include noise mitigation 58-60 , quantum random access memory [61][62][63] , tensor networks 64 , Ansatz architecture, non-linear algorithms and cost-efficient encoding.

Data availability
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

A Fully implicit scheme for linear Hermitian source term
For a reaction-diffusion system of equations with non-trivial source terms (35), time-stepping could be rendered fully implicit if the source terms can be expressed as a linear Hermitian matrix with constant coefficients. Consider a two-component 1D chemical reaction with source terms of the form where k ij ( i, j ∈ {1, 2} ) are kinetic rate constants that are components of the linear Hermitian matrix K = K T ∈ R N×N , whose off-diagonal elements are equal, i.e. k 12 = k 21 ; I is the identity matrix of size N × N and [u 1 , u 2 ] T is a concentration vector of length 2N. This problem requires n + 1 qubits, where n = log 2 N . Following the implicit Euler scheme (Eq. 8), we set up a 2N × 2N coefficient matrix, which decomposes as where A x,β is the discretized N × N coefficient matrix for a single component (Eq. 16) containing up to four Hamiltonian terms H 1−4 . Note the last three additional Hamiltonian terms contributed by the source term.

B Vorticity stream-function formulation
For two-dimensional incompressible flows, the vorticity stream-function formulation can be used to eliminate the pressure as a dependent variable, such that where ω = ∂v/∂x − ∂u/∂y is the flow vorticity and the stream-function ψ satisfies u = −∂ψ/∂y and v = ∂ψ/∂x . It follows that the re-formulation {u, v, p} → {ω, ψ} can simplify the variational quantum algorithm 3 by reducing the number of velocity components by one, at the cost of specifying stream-function values along the domain boundaries. (A.2) A = I ⊗n+1 + δ x I ⊗ A x,β − �t(k 11 I 0 + k 22 I 1 + k 12 X) ⊗ I ⊗n ,